Take-home Exercise 03

Author

Andrea Yeo

Published

March 22, 2025

Modified

March 27, 2025

Prototyping Modules for Visual Analytics Shiny Application

1. The task

For this take-home exercise, each team member have to choose one module from our proposed Shiny application and complete the following tasks:

  1. Verify package support: Identify and ensure all required R packages for the module are available on CRAN.
  2. Code testing: Prepare and test the relevant R code to confirm it runs correctly and produces the expected output.
  3. Define inputs and outputs: Specify the parameters (inputs) and outputs that will be featured in the Shiny application.
  4. Select UI components: Choose suitable Shiny UI elements to present and control the identified parameters within the app interface.

2. Getting started

Our project utilizes open-source data from the World Happiness Report – Data Sharing covering the years 2015 to 2025. Detailed information on the data cleaning and preparation process is available on our group’s Netlify site, under the “Data Preparation” tab: Group Visual Analytic Application.

Specifically, our project focuses on the visualisation of global happiness trends. I will be exploring the modules related to geospatial and aspatial analysis to explore both spatial patterns and non-spatial factors influencing happiness scores across countries. For the prototyping of these modules, I will be using the cleaned_data dataset prepared after the data processing phase.

2.1 Loading R packages

The below R packages will be used in this exercise and for the Shiny application

Code 
pacman::p_load(sf, tidyverse, tmap, dplyr,
               spatstat, spdep,
               lubridate, leaflet,
               plotly, DT, viridis,
               ggplot2, sfdep, janitor, shiny, rnaturalearth)

2.2 Reading the dataset

Code 
happiness <- read.csv("data/world_happiness.csv")

2.3 Understanding the data structure

  • glimpse(): provides a transposed overview of a dataset, showing variables and their types in a concise format.
  • head(): displays the first few rows of a dataset (default is 6 rows) to give a quick preview of the data.
  • summary(): generates a statistical summary of each variable, including measures like mean, median, and range for numeric data.
  • duplicated():returns a logical vector indicating which elements or rows in a vector or data frame are duplicates.
  • colSums(is.na()): counts the number of missing values (NA) in each column of the data frame.
  • str(): use str() to display the column names, data types, and a preview of the data.
Code 
glimpse(happiness)
Rows: 1,656
Columns: 12
$ year                 <int> 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2…
$ country              <chr> "Afghanistan", "Afghanistan", "Afghanistan", "Afg…
$ ladder_score         <dbl> 3.5750, 3.3600, 3.7940, 3.6320, 3.2030, 2.5669, 2…
$ economy_score        <dbl> 0.31982, 0.38227, 0.40148, 0.33200, 0.35000, 0.30…
$ social_score         <dbl> 0.30285, 0.11037, 0.58154, 0.53700, 0.51700, 0.35…
$ lifeexpectancy_score <dbl> 0.30335, 0.17344, 0.18075, 0.25500, 0.36100, 0.26…
$ freedom_score        <dbl> 0.23414, 0.16430, 0.10618, 0.08500, 0.00000, 0.00…
$ generosity_score     <dbl> 0.36510, 0.31268, 0.31187, 0.19100, 0.15800, 0.13…
$ corrperception_score <dbl> 0.09719, 0.07112, 0.06116, 0.03600, 0.02500, 0.00…
$ residual_score       <dbl> 1.95210, 2.14558, 2.15080, 2.19600, 1.79200, 1.50…
$ rank                 <int> 153, 154, 141, 145, 154, 153, 149, 146, 137, 143,…
$ region               <chr> "Asia", "Asia", "Asia", "Asia", "Asia", "Asia", "…
Code 
head(happiness)
  year     country ladder_score economy_score social_score lifeexpectancy_score
1 2014 Afghanistan       3.5750       0.31982      0.30285              0.30335
2 2015 Afghanistan       3.3600       0.38227      0.11037              0.17344
3 2016 Afghanistan       3.7940       0.40148      0.58154              0.18075
4 2017 Afghanistan       3.6320       0.33200      0.53700              0.25500
5 2018 Afghanistan       3.2030       0.35000      0.51700              0.36100
6 2019 Afghanistan       2.5669       0.30071      0.35643              0.26605
  freedom_score generosity_score corrperception_score residual_score rank
1       0.23414          0.36510              0.09719        1.95210  153
2       0.16430          0.31268              0.07112        2.14558  154
3       0.10618          0.31187              0.06116        2.15080  141
4       0.08500          0.19100              0.03600        2.19600  145
5       0.00000          0.15800              0.02500        1.79200  154
6       0.00000          0.13523              0.00123        1.50724  153
  region
1   Asia
2   Asia
3   Asia
4   Asia
5   Asia
6   Asia
Code 
summary(happiness)
      year        country           ladder_score   economy_score   
 Min.   :2014   Length:1656        Min.   :1.364   Min.   :0.0000  
 1st Qu.:2016   Class :character   1st Qu.:4.606   1st Qu.:0.7607  
 Median :2019   Mode  :character   Median :5.497   Median :1.0983  
 Mean   :2019                      Mean   :5.460   Mean   :1.0763  
 3rd Qu.:2022                      3rd Qu.:6.303   3rd Qu.:1.3946  
 Max.   :2024                      Max.   :7.842   Max.   :2.2090  
  social_score    lifeexpectancy_score freedom_score    generosity_score
 Min.   :0.0000   Min.   :0.0000       Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.8642   1st Qu.:0.4028       1st Qu.:0.3794   1st Qu.:0.1069  
 Median :1.1122   Median :0.5963       Median :0.4972   Median :0.1700  
 Mean   :1.0782   Mean   :0.5756       Mean   :0.4914   Mean   :0.1846  
 3rd Qu.:1.3496   3rd Qu.:0.7599       3rd Qu.:0.6110   3rd Qu.:0.2421  
 Max.   :1.8400   Max.   :1.1410       Max.   :1.0180   Max.   :0.8381  
 corrperception_score residual_score        rank           region         
 Min.   :0.0000       Min.   :-0.110   Min.   :  1.00   Length:1656       
 1st Qu.:0.0590       1st Qu.: 1.553   1st Qu.: 38.00   Class :character  
 Median :0.1010       Median : 1.922   Median : 76.00   Mode  :character  
 Mean   :0.1356       Mean   : 1.919   Mean   : 75.97                     
 3rd Qu.:0.1726       3rd Qu.: 2.309   3rd Qu.:113.25                     
 Max.   :0.5870       Max.   : 3.838   Max.   :158.00
Code 
happiness[duplicated(happiness),]
 [1] year                 country              ladder_score        
 [4] economy_score        social_score         lifeexpectancy_score
 [7] freedom_score        generosity_score     corrperception_score
[10] residual_score       rank                 region              
<0 rows> (or 0-length row.names)
Note
  • Ensure that there are no duplicated columns, if not will have to investigate further.
Code 
colSums(is.na(happiness))
                year              country         ladder_score 
                   0                    0                    0 
       economy_score         social_score lifeexpectancy_score 
                   0                    0                    0 
       freedom_score     generosity_score corrperception_score 
                   0                    0                    0 
      residual_score                 rank               region 
                   0                    0                    0
Note
  • Ensure that there are no NA values, if not will have to investigate further.
  • Possibility to use drop_na() function to drop rows where any specified column contains a missing value.
Code 
str(happiness)
'data.frame':   1656 obs. of  12 variables:
 $ year                : int  2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 ...
 $ country             : chr  "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
 $ ladder_score        : num  3.58 3.36 3.79 3.63 3.2 ...
 $ economy_score       : num  0.32 0.382 0.401 0.332 0.35 ...
 $ social_score        : num  0.303 0.11 0.582 0.537 0.517 ...
 $ lifeexpectancy_score: num  0.303 0.173 0.181 0.255 0.361 ...
 $ freedom_score       : num  0.234 0.164 0.106 0.085 0 ...
 $ generosity_score    : num  0.365 0.313 0.312 0.191 0.158 ...
 $ corrperception_score: num  0.0972 0.0711 0.0612 0.036 0.025 ...
 $ residual_score      : num  1.95 2.15 2.15 2.2 1.79 ...
 $ rank                : int  153 154 141 145 154 153 149 146 137 143 ...
 $ region              : chr  "Asia" "Asia" "Asia" "Asia" ...
Note
  • Ensure that all variables are correctly classified by data type; recast variable types if needed.
  • Variables are correctly classified - where categorical variables are classified as character, while continuous variables are classified as double.

The happiness tibble contains 12 attributes, as shown above.

The following preprocessing checks were conducted as part of data preparation:

Preprocessing Checks
  • Verified that the correct data types were loaded in the happiness dataset using glimpse() and str()
  • Ensured there were no duplicate variable names using duplicated() in the dataset
  • Checked for missing values using colSums(is.na())
Data scale review
  • The column values are already normalized, ranging consistently between 1 and 10 across countries.
  • No further standardization/ transformation is required, as the scale is uniform and interpretable for spatial analysis.

3. Geospatial analysis (Choropleth Map)

3.1 To load world country boundaries

To support the geospatial visualisation of global happiness trends, the rnaturalearth package was used to import spatial boundary data of countries. Specifically, the ne_countries() function was called with a medium scale and returned as a simple features (sf) object, stored in the world variable. This spatial dataset provides the geographic outlines of all countries and is essential for merging with the World Happiness Report data to enable choropleth mapping and spatial analysis at the country level.

Code 
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

3.2 Join happiness data with geospatial features

To prepare the dataset for mapping, the World Happiness data was filtered to include only records from the year 2024. This filtered data (happiness_latest) was then joined with the spatial dataset (world) using the country name as the common key. The resulting dataset, world_happy, combines both geographic boundaries and happiness scores, enabling geospatial visualisation of happiness levels across countries.

Code 
# Join by country name 
happiness_latest <- happiness %>% 
  filter(year == 2024)  

world_happy <- world %>%
  left_join(happiness_latest, by = c("name" = "country"))

3.3 Choropleth map of global happiness (basic version)

The code below plots a choropleth map to visualise the distribution of happiness scores across countries for the year 2024. Using the tmap package in plotting mode, the tm_shape() function defines the spatial data object (world_happy), which contains both geographic boundaries and the associated happiness data. The tm_polygons() function maps the ladder_score variable to a blue color gradient, with darker shades indicating higher happiness scores.

Finally, tm_layout() adds a descriptive title and positions the legend outside the map area to improve clarity and presentation.

Code 
tmap_mode("plot")
tmap_style("white")

tm_shape(world_happy) +
  tm_polygons("ladder_score", palette = "Blues", title = "Happiness Score (2024)") +
  tm_layout(title = "World Happiness Map", legend.outside = TRUE)

3.4 Choropleth map with styling and classification (enhanced version)

The enhanced choropleth map uses the tmap package with improved styling and classification. It applies a quantile-based scheme with five classes and the "YlOrRd" palette to represent ladder_score, where darker shades indicate higher happiness. Additional cartographic elements such as a compass, scale bar, borders, and descriptive credits are included to improve readability and interpretability. The legend is placed on the right, and a classic map style is applied for a cleaner layout.

Code 
tmap_style("classic")
tmap_mode("plot")

tm_shape(world_happy) +
  tm_polygons(
    col = "ladder_score",
    palette = "YlOrRd",
    style = "quantile",
    n = 5,
    title = "Happiness Score (2024)"
  ) +
  tm_layout(
    main.title = "World Happiness Distribution (2024)",
    main.title.position = "center",
    main.title.size = 1.2,
    legend.outside = TRUE,
    legend.position = c("right", "center"),
    frame = FALSE
  ) +
  tm_borders(alpha = 0.4, lwd = 0.3) +
  tm_compass(type = "8star", size = 2, position = c("left", "top")) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_credits("Source: World Happiness Report 2024", position = c("left", "bottom")) +
  tm_credits("Note: Darker shades indicate higher happiness; lighter shades indicate lower happiness.",
             position = c("left", "bottom"), just = "right", size = 0.6)

3.5 Summary of differences between Basic version and Enhanced Version

Feature Basic Version Enhanced Version
Classification method Default(“pretty”) Custom(“quantile”, 5 bins)
Layout Basic title, default layout Custom title, centered, larger size
Map Borders None Added with transparency and thin width
Compass/ Scale Bar Not included Included
Credits/ Source Not included Included
Frame Default Removed

4. Drawing small multiple (faceted) map

The code below creates a faceted choropleth map to visualise global happiness scores from 2014 to 2024. Using tmap, each facet represents a different year, allowing for easy comparison of happiness trends across countries. Geospatial data is joined using rnaturalearth.

Code 
# Load required libraries
library(tidyverse); library(tmap); library(sf); library(rnaturalearth)

# Set style and plotting mode
tmap_style("white")
tmap_mode("plot")

# Step 1: Filter for all years from 2014 to 2024
years_to_plot <- 2014:2024

happiness_subset <- happiness %>%
  filter(year %in% years_to_plot)

# Step 2: Load and join with world spatial data
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

world_happy_facet <- world %>%
  left_join(happiness_subset, by = c("name" = "country"))

# Step 3: Filter out rows with missing ladder_score
world_happy_facet_clean <- world_happy_facet %>%
  filter(!is.na(ladder_score))

# Step 4: Plot small multiple choropleth maps
tm_shape(world_happy_facet_clean) +
  tm_polygons(
    col = "ladder_score",
    palette = "YlGnBu",
    style = "quantile",
    title = "Happiness Score"
  ) +
  tm_facets(by = "year", ncol = 2) +
  tm_layout(
    main.title = "Global Happiness Scores (2014–2024)",
    main.title.position = "center",
    panel.label.size = 2.0,
    legend.outside = TRUE,
    legend.position = c("center", "top"),   
    legend.height = 1.5,
    legend.width = 1.0,
    legend.text.size = 1.5,
    legend.title.size = 1.5,
    asp = 0,
    between.margin = 0,  
    outer.margins = c(0, 0, 0, 0)  
  )

5. Animating global happiness over time

The code below produces an animated choropleth map using ggplot2, gganimate, and sf to illustrate the temporal evolution of global happiness scores from 2014 to 2024. Country-level geospatial data from rnaturalearth is merged with World Happiness Report data via a spatial join. transition_time() is used to animate by year, with countries shaded according to their ladder score.

Note

Note: Although the underlying world geometries are static, the animated map may appear to shift or jitter between frames. This occurs because some countries are missing data in certain years, leading to variations in the number of plotted geometries. As a result, gganimate recalculates the plot layout for each frame, even when coord_sf(xlim, ylim, expand = FALSE) is specified. The visual movement is not due to changes in the map itself, but rather to inconsistencies in data availability across years.

Code 
library(tidyverse); library(sf); library(gganimate); library(rnaturalearth)

# Load base world map
world <- ne_countries(scale = "medium", returnclass = "sf") %>%
  select(name, geometry)

# Load and clean happiness data
happiness <- read_csv("data/world_happiness.csv") %>%
  mutate(year = as.numeric(year)) %>%
  filter(year %in% 2014:2024)

# Merge first — keep all geometries
world_happy <- world %>%
  left_join(happiness, by = c("name" = "country")) %>%
  mutate(year = as.numeric(year))  

# Plot
ggplot(world_happy) +
  geom_sf(aes(fill = ladder_score)) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1, na.value = "lightgray") +
  labs(title = "Happiness Score: {frame_time}", fill = "Score") +
  transition_time(year) +
  ease_aes('linear') +
  coord_sf(
    xlim = c(-180, 180),
    ylim = c(-60, 90),
    expand = FALSE
  )

6. Exploring global happiness with an interactive map - Choropleth Map

This code creates an interactive world map using tmap, showing happiness scores by country for 2024. Each country is colored based on its ladder score, and tooltips display the country name when hovered. The map allows users to zoom and explore happiness data visually.

Code 
happiness_latest <- happiness_subset %>%
  filter(year == 2024)

world_happy_latest <- world %>%
  left_join(happiness_latest, by = c("name" = "country"))

world_happy_latest_clean <- world_happy_latest %>%
  filter(!is.na(ladder_score))

tmap_mode("view")

tm_shape(world_happy_latest_clean) +
  tm_polygons(
    col = "ladder_score",
    palette = "YlGnBu",
    id = "name",  # this is what shows in the tooltip
    title = "Happiness Score (2024)"
  )

7. User interface design - Part 1 - Choropleth Map

In this module, I developed an interactive Shiny application to visualize the World Happiness Score by country and year. The key input is a dropdown menu (selectInput) that allows users to choose a specific year from 2014 to 2024. Based on the selected year, the output is a dynamic choropleth map (tmapOutput) showing each country’s happiness score using a color gradient.

Note

  • Input: A single selectInput() allows the user to choose the Year of interest from 2014 to 2024.

  • Output: A dynamic choropleth map is rendered using tmap, updating automatically based on the selected year.

  • UI Components Used: selectInput() for the dropdown menu, sidebarPanel() for input placement, and tmapOutput() in mainPanel() for rendering the interactive map.

Code 
# Load libraries
library(shiny); library(tmap); library(tidyverse); library(sf); library(rnaturalearth)

# Prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
  mutate(year = as.numeric(year))

years_to_plot <- 2014:2024
happiness_subset <- happiness %>%
  filter(year %in% years_to_plot)

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

world_happy <- world %>%
  left_join(happiness_subset, by = c("name" = "country"))

# UI
ui <- fluidPage(
  titlePanel("🌍 World Happiness Score by Year"),
  
  sidebarLayout(
    sidebarPanel(
      selectInput("selected_year", "Select Year:",
                  choices = sort(unique(world_happy$year)), 
                  selected = 2024)
    ),
    
    mainPanel(
      tmapOutput("happiness_map")
    )
  )
)

# Server
server <- function(input, output, session) {
  
  filtered_data <- reactive({
    world_happy %>%
      filter(year == input$selected_year & !is.na(ladder_score))
  })
  
  output$happiness_map <- renderTmap({
    tmap_mode("view")
    tm_shape(filtered_data()) +
      tm_polygons(
        col = "ladder_score",
        palette = "YlGnBu",
        id = "name",
        popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
        title = paste("Happiness Score:", input$selected_year)
      )
  })
}

# Run the app
shinyApp(ui, server)

Shiny applications not supported in static R Markdown documents

8. Mapping global happiness with proportional circles - Proportional Symbol Map

This code below generates an interactive map using leaflet to visualize 2024 World Happiness Scores. Each country is displayed as a circle located at its geographic center, with the size representing its happiness score. I added tooltips which provide additional insights like economy, life expectancy, freedom, and region, and a legend which helps interpret the happiness scores visually.

Code 
library(tidyverse); library(readr); library(leaflet)
library(sf); library(rnaturalearth); library(rnaturalearthdata)

# Load the world happiness data
happiness <- read_csv("data/world_happiness.csv")

# Filter for a specific year (e.g., 2024)
happiness_latest <- happiness %>%
  filter(year == 2024)

# Load country polygons and calculate centroids
world_sf <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
  select(name, geometry)

# Join happiness data to spatial data
world_happy_sf <- world_sf %>%
  left_join(happiness_latest, by = c("name" = "country")) %>%
  filter(!is.na(ladder_score))

# Get centroids for each country
world_happy_points <- st_centroid(world_happy_sf)

# Extract coordinates for leaflet
world_happy_points_coords <- world_happy_points %>%
  mutate(
    lon = st_coordinates(geometry)[, 1],
    lat = st_coordinates(geometry)[, 2]
  )

leaflet(world_happy_points_coords) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addCircleMarkers(
    lng = ~lon,
    lat = ~lat,
    radius = ~ladder_score * 3,
    color = "black",
    fillColor = ~colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score)(ladder_score),
    fillOpacity = 0.6,  
    stroke = TRUE,
    weight = 0.5,
    popup = ~paste0(
      "<b>Country:</b> ", name, "<br/>",
      "<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
      "<b>Economy:</b> ", round(economy_score, 2), "<br/>",
      "<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
      "<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
      "<b>Region:</b> ", region
    )
  ) %>%
  addLegend(
    "bottomright",
    pal = colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score),
    values = ~ladder_score,
    title = "Happiness Score",
    opacity = 1
  )

9. User interface design - Part 2 - Proportional Symbol Map

In this module, I developed an interactive Shiny application that allows users to explore the World Happiness Score by country for the year 2024. The main input is a country selector using selectInput, which enables users to search and select a country from the dropdown list. Once selected, the output is an interactive leaflet map that automatically zooms to the selected country and displays a popup with additional details such as economy, life expectancy, freedom, and region.

Note

  • Input: A selectInput() enables users to search and select a country from the list of available options for 2024.

  • Output: A proportional symbol map rendered via leaflet, displaying country-level happiness scores using interactive circle markers with popups.

  • UI Components Used: selectInput() in the sidebarPanel() for country selection, and leafletOutput() in the mainPanel() for displaying the dynamic map.

Code 
library(shiny); library(tidyverse); library(readr); library(leaflet)
library(sf); library(rnaturalearth); library(rnaturalearthdata)

# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
  filter(year == 2024)

world_sf <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
  select(name, geometry)

world_happy_sf <- world_sf %>%
  left_join(happiness, by = c("name" = "country")) %>%
  filter(!is.na(ladder_score))

world_happy_points <- st_centroid(world_happy_sf)

world_happy_points_coords <- world_happy_points %>%
  mutate(
    lon = st_coordinates(geometry)[, 1],
    lat = st_coordinates(geometry)[, 2]
  )

# UI
ui <- fluidPage(
  titlePanel("World Happiness (2024)"),
  sidebarLayout(
    sidebarPanel(
      selectInput("selected_country", "Search Country:", 
                  choices = sort(world_happy_points_coords$name),
                  selected = NULL)
    ),
    mainPanel(
      leafletOutput("map", height = "400px")
    )
  )
)

# Server
server <- function(input, output, session) {
  
  output$map <- renderLeaflet({
    leaflet(world_happy_points_coords) %>%
      addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
      addCircleMarkers(
        lng = ~lon,
        lat = ~lat,
        radius = ~ladder_score * 3,
        color = "black",
        fillColor = ~colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score)(ladder_score),
        fillOpacity = 0.6,
        stroke = TRUE,
        weight = 0.5,
        popup = ~paste0(
          "<b>Country:</b> ", name, "<br/>",
          "<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
          "<b>Economy:</b> ", round(economy_score, 2), "<br/>",
          "<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
          "<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
          "<b>Region:</b> ", region
        )
      ) %>%
      addLegend(
        "bottomright",
        pal = colorNumeric("YlGnBu", domain = world_happy_points_coords$ladder_score),
        values = ~ladder_score,
        title = "Happiness Score",
        opacity = 1
      )
  })
  
  observeEvent(input$selected_country, {
    selected_data <- world_happy_points_coords %>%
      filter(name == input$selected_country)
    
    leafletProxy("map") %>%
      setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
      clearPopups() %>%
      addPopups(
        lng = selected_data$lon,
        lat = selected_data$lat,
        popup = paste0(
          "<b>Country:</b> ", selected_data$name, "<br/>",
          "<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
          "<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
          "<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
          "<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
          "<b>Region:</b> ", selected_data$region
        )
      )
  })
}

# Run the app
shinyApp(ui, server)

Shiny applications not supported in static R Markdown documents

10. User interface design - Part 3 - Combining both Choropleth Map & Proportional Symbol Map

This module presents an interactive Shiny application with two coordinated maps to explore World Happiness data. The Choropleth Map allows users to select a year and view global happiness patterns using a color gradient. In contrast, the Proportional Symbol Map lets users search for a country, zoom in automatically, and view detailed indicators such as economy, life expectancy, and freedom. Circle size and color reflect the happiness score.

After considering various approaches, I chose to implement both map types as they serve different purposes:

  • the choropleth reveals macro-level trends, while

  • the proportional map provides micro-level insights.

Tip
  • Note: This application includes smart interactivity features. Selecting a country from the dropdown automatically zooms both maps to the selected country for focused viewing. Conversely, clicking a country on the choropleth map will update the dropdown menu —ensuring two-way interaction. These features improve usability by making geographic exploration more intuitive and dynamic.
Note

  • Input: Two selectInput() components allow users to:

    • (1) choose the year (2014–2024), and

    • (2) search and select a country for focused exploration.

  • Output: A coordinated view using two maps:

    • A choropleth map rendered via tmap, showing happiness levels by country with zoom-to-country feature.

    • A proportional symbol map rendered with leaflet, displaying happiness scores with interactive popups.

  • UI Components Used: selectInput() (for year and country), sidebarPanel() for inputs, and tmapOutput() + leafletOutput() within a fluidRow() layout for side-by-side map rendering.

Code 
# Load libraries
library(shiny); library(tidyverse); library(sf); library(tmap); library(leaflet); library(rnaturalearth)

# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
  mutate(year = as.numeric(year)) %>%
  filter(year %in% 2014:2024)

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

world_happy <- world %>%
  left_join(happiness, by = c("name" = "country"))

ui <- fluidPage(
  titlePanel("🌍 World Happiness Explorer"),
  sidebarLayout(
    sidebarPanel(
      selectInput("selected_year", "Select Year:", choices = sort(unique(world_happy$year)), selected = 2024),
      selectInput("selected_country", "Search Country:", choices = NULL)
    ),
    mainPanel(
      fluidRow(
        column(6,
               h4("Choropleth Map"),
               tmapOutput("choropleth_map", height = "500px")
        ),
        column(6,
               h4("Proportional Symbol Map"),
               leafletOutput("prop_map", height = "500px")
        )
      )
    )
  )
)

server <- function(input, output, session) {

  # Reactive filtered data by year
  filtered_data <- reactive({
    world_happy %>%
      filter(year == input$selected_year & !is.na(ladder_score))
  })

  # Update dropdown
  observe({
    updateSelectInput(session, "selected_country",
                      choices = sort(unique(filtered_data()$name)))
  })

  # Choropleth Map - Zoom to selected country
  output$choropleth_map <- renderTmap({
    tmap_mode("view")

    # Zooming to selected country
    selected_geom <- filtered_data() %>% filter(name == input$selected_country)
    bbox_zoom <- if (nrow(selected_geom) > 0) st_bbox(selected_geom) else st_bbox(filtered_data())

    tm_shape(filtered_data(), bbox = bbox_zoom) +
      tm_polygons(
        col = "ladder_score",
        palette = "YlGnBu",
        id = "name",
        popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
        title = paste("Happiness Score:", input$selected_year)
      )
  })

  # Proportional Symbol Map
  output$prop_map <- renderLeaflet({
    centroids <- st_centroid(filtered_data())
    coords <- centroids %>%
      mutate(
        lon = st_coordinates(geometry)[, 1],
        lat = st_coordinates(geometry)[, 2]
      )

    leaflet(coords) %>%
      addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
      addCircleMarkers(
        lng = ~lon,
        lat = ~lat,
        radius = ~ladder_score * 3,
        color = "black",
        fillColor = ~colorNumeric("YlGnBu", domain = coords$ladder_score)(ladder_score),
        fillOpacity = 0.6,
        stroke = TRUE,
        weight = 0.5,
        popup = ~paste0(
          "<b>Country:</b> ", name, "<br/>",
          "<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
          "<b>Economy:</b> ", round(economy_score, 2), "<br/>",
          "<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
          "<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
          "<b>Region:</b> ", region
        )
      ) %>%
      addLegend(
        "bottomright",
        pal = colorNumeric("YlGnBu", domain = coords$ladder_score),
        values = ~ladder_score,
        title = "Happiness Score",
        opacity = 1
      )
  })

  # Auto-zoom on dropdown change
  observeEvent(input$selected_country, {
    centroids <- st_centroid(filtered_data())
    coords <- centroids %>%
      mutate(
        lon = st_coordinates(geometry)[, 1],
        lat = st_coordinates(geometry)[, 2]
      )

    selected_data <- coords %>% filter(name == input$selected_country)

    leafletProxy("prop_map") %>%
      setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
      clearPopups() %>%
      addPopups(
        lng = selected_data$lon,
        lat = selected_data$lat,
        popup = paste0(
          "<b>Country:</b> ", selected_data$name, "<br/>",
          "<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
          "<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
          "<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
          "<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
          "<b>Region:</b> ", selected_data$region
        )
      )
  })

  # Auto-update dropdown when user clicks map
  observeEvent(input$choropleth_map_shape_click, {
    clicked_country <- input$choropleth_map_shape_click$id
    updateSelectInput(session, "selected_country", selected = clicked_country)
  })
}

shinyApp(ui, server)

Shiny applications not supported in static R Markdown documents

11. User interface design - Part 4 - Adding regions

In this iteration, I introduced an additional selectInput() to filter the data by region. This allows users to explore happiness scores by broader geographical areas before narrowing down to a specific country. The country dropdown is dynamically updated based on the selected region, and selecting a country will also auto-update the region field—ensuring a smooth and linked filtering experience.

Tip
  • Note: To explore and ensure that the region and country drop downs are dynamically linked, where
    • selecting a region will show only countries in that region
    • selecting a country will auto-update the region selection to match
Note

  • Input: Three selectInput() components allow users to:

    • (1) choose the year (2014–2024),

    • (2) filter by region,

    • (3) search and select a country within a selected region.

  • Output: A coordinated view using two maps:

    • A choropleth map rendered via tmap, showing happiness levels by country with zoom-to-country feature.

    • A proportional symbol map rendered with leaflet, displaying happiness scores with interactive popups.

  • UI Components Used: selectInput() (for year, region and country), sidebarPanel() for inputs, and tmapOutput() + leafletOutput() within a fluidRow() layout for side-by-side map rendering.

Code 
# Load libraries
library(shiny); library(tidyverse); library(sf); library(tmap); library(leaflet); library(rnaturalearth); library(readr)

# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
  mutate(year = as.numeric(year)) %>%
  filter(year %in% 2014:2024)

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world_happy <- world %>%
  left_join(happiness, by = c("name" = "country"))

# UI
ui <- fluidPage(
  titlePanel("🌍 World Happiness Explorer"),
  sidebarLayout(
    sidebarPanel(
      selectInput("selected_year", "Select Year:",
                  choices = sort(unique(world_happy$year)),
                  selected = 2024),
      selectInput("selected_region", "Filter by Region:", choices = NULL),
      selectInput("selected_country", "Search Country:", choices = NULL)
    ),
    mainPanel(
      fluidRow(
        column(6,
               h4("Choropleth Map"),
               tmapOutput("choropleth_map", height = "500px")
        ),
        column(6,
               h4("Proportional Symbol Map"),
               leafletOutput("prop_map", height = "500px")
        )
      )
    )
  )
)

# Server
server <- function(input, output, session) {

  # Initial values for dropdowns
  observe({
    regions <- sort(unique(world_happy$region))
    updateSelectInput(session, "selected_region",
                      choices = c("All", regions),
                      selected = "All")
  })

  # Reactive filtered data
  filtered_data <- reactive({
    data <- world_happy %>%
      filter(year == input$selected_year & !is.na(ladder_score))
    if (input$selected_region != "All") {
      data <- data %>% filter(region == input$selected_region)
    }
    return(data)
  })

  # Update country dropdown when region changes
  observeEvent(input$selected_region, {
    countries <- filtered_data() %>% pull(name) %>% unique() %>% sort()
    updateSelectInput(session, "selected_country",
                      choices = countries,
                      selected = countries[1])
  })

  # Update region when country changes
  observeEvent(input$selected_country, {
    selected_region <- world_happy %>%
      filter(name == input$selected_country,
             year == input$selected_year) %>%
      pull(region) %>% unique()

    if (!is.null(selected_region) && length(selected_region) == 1) {
      updateSelectInput(session, "selected_region",
                        selected = selected_region)
    }
  })

  # Choropleth Map
  output$choropleth_map <- renderTmap({
    tmap_mode("view")
    selected_geom <- filtered_data() %>% filter(name == input$selected_country)
    bbox_zoom <- if (nrow(selected_geom) > 0) st_bbox(selected_geom) else st_bbox(filtered_data())

    tm_shape(filtered_data(), bbox = bbox_zoom) +
      tm_polygons(
        col = "ladder_score",
        palette = "YlGnBu",
        id = "name",
        popup.vars = c("Country" = "name", "Happiness" = "ladder_score"),
        title = paste("Happiness Score:", input$selected_year)
      )
  })

  # Proportional Symbol Map
  output$prop_map <- renderLeaflet({
    centroids <- st_centroid(filtered_data())
    coords <- centroids %>%
      mutate(
        lon = st_coordinates(geometry)[, 1],
        lat = st_coordinates(geometry)[, 2]
      )

    leaflet(coords) %>%
      addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
      addCircleMarkers(
        lng = ~lon,
        lat = ~lat,
        radius = ~ladder_score * 3,
        color = "black",
        fillColor = ~colorNumeric("YlGnBu", domain = coords$ladder_score)(ladder_score),
        fillOpacity = 0.6,
        stroke = TRUE,
        weight = 0.5,
        popup = ~paste0(
          "<b>Country:</b> ", name, "<br/>",
          "<b>Happiness Score:</b> ", round(ladder_score, 2), "<br/>",
          "<b>Economy:</b> ", round(economy_score, 2), "<br/>",
          "<b>Life Expectancy:</b> ", round(lifeexpectancy_score, 2), "<br/>",
          "<b>Freedom:</b> ", round(freedom_score, 2), "<br/>",
          "<b>Region:</b> ", region
        )
      ) %>%
      addLegend(
        "bottomright",
        pal = colorNumeric("YlGnBu", domain = coords$ladder_score),
        values = ~ladder_score,
        title = "Happiness Score",
        opacity = 1
      )
  })

  # Auto-zoom and popup on country selection
  observeEvent(input$selected_country, {
    centroids <- st_centroid(filtered_data())
    coords <- centroids %>%
      mutate(
        lon = st_coordinates(geometry)[, 1],
        lat = st_coordinates(geometry)[, 2]
      )

    selected_data <- coords %>% filter(name == input$selected_country)

    if (nrow(selected_data) > 0 &&
        !is.na(selected_data$lon) &&
        !is.na(selected_data$lat)) {
      leafletProxy("prop_map") %>%
        setView(lng = selected_data$lon, lat = selected_data$lat, zoom = 5) %>%
        clearPopups() %>%
        addPopups(
          lng = selected_data$lon,
          lat = selected_data$lat,
          popup = paste0(
            "<b>Country:</b> ", selected_data$name, "<br/>",
            "<b>Happiness Score:</b> ", round(selected_data$ladder_score, 2), "<br/>",
            "<b>Economy:</b> ", round(selected_data$economy_score, 2), "<br/>",
            "<b>Life Expectancy:</b> ", round(selected_data$lifeexpectancy_score, 2), "<br/>",
            "<b>Freedom:</b> ", round(selected_data$freedom_score, 2), "<br/>",
            "<b>Region:</b> ", selected_data$region
          )
        )
    }
  })
}

# Run the app
shinyApp(ui, server)

Shiny applications not supported in static R Markdown documents

12. Exploring the LISA cluster map

In the next tab of our Cluster and Outlier Analysis module, I will apply spatial autocorrelation techniques to the World Happiness dataset by generating a Moran Scatter Plot and a Local Indicators of Spatial Association (LISA) Cluster Map.

The Moran Scatter Plot will quantify global spatial autocorrelation, indicating the overall degree of spatial clustering of happiness scores. The LISA Cluster Map will identify statistically significant local clusters—such as High-High (hot spots), Low-Low (cold spots), and spatial outliers (High-Low, Low-High)—allowing for granular analysis of localized spatial dependencies. T

12.1 Plotting the Moran Scatter Plot

The Moran scatter plot visualizes the spatial autocorrelation between countries’ Happiness Scores and the average scores of their neighboring countries. This helps us understand whether similar levels of happiness cluster geographically.

In the example below, we use moran.plot() from the spdep package to explore spatial relationships in the 2024 Happiness Score:

Interpretation:

The plot is split into four quadrants, each representing a different spatial relationship:

  • High-High (HH) - Top right quadrant: countries with high

  • Low-Low (LL) - Bottom-left quadrant: countries with low happiness surrounded by similarly unhappy neighbors.

  • High-Low (HL) - Bottom-right quadrant: — countries with high happiness surrounded by low-scoring neighbors (potential outliers).

  • Low-High (LH) - Top-left: countries with low happiness surrounded by happy neighbors (also outliers).

Tip

Note: Scatter plot provides an intuitive way to detect potential spatial clusters/ anomalies.

Code 
# Filter and prepare data
world_2024 <- world_happy %>% filter(year == 2024 & !is.na(ladder_score))

# Create spatial weights (Queen's case)
coords <- st_coordinates(st_centroid(world_2024))
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Moran Scatter Plot
moran.plot(world_2024$ladder_score, lw,
           labels = as.character(world_2024$name),
           xlab = "Happiness Score (2024)",
           ylab = "Spatially Lagged Happiness Score",
           zero.policy = TRUE)

12.2 Preparing LISA cluster class for Word Happiness (2024)

Before generating the LISA cluster map, it is important to center both the spatially lagged variable and the local Moran’s I values. This step allows us to identify meaningful spatial patterns by determining whether a country’s happiness score — and the scores of its neighbors — are above or below the global mean.

By doing so, we can categorize each country into one of the four LISA cluster types (High-High, Low-Low, High-Low, Low-High). These classifications are only meaningful when assessed relative to the dataset’s overall distribution, making centering a crucial step for accurate spatial interpretation.

Code 
# Load required library
library(spdep)

# Step 1: Filter 2024 data
world_2024 <- world_happy %>% 
  filter(year == 2024 & !is.na(ladder_score))

# Step 2: Create spatial weights
coords <- st_coordinates(st_centroid(world_2024))
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Step 3: Calculate local Moran's I
localMI <- localmoran(world_2024$ladder_score, lw, zero.policy = TRUE)

# Step 4: Spatial lag of the variable
world_2024$lag_score <- lag.listw(lw, world_2024$ladder_score)

# Step 5: Center spatial lag & local Moran's I
DV <- world_2024$lag_score - mean(world_2024$lag_score, na.rm = TRUE)
LM_I <- localMI[, 1] - mean(localMI[, 1], na.rm = TRUE)

# Step 6: Set significance threshold
signif <- 0.05

# Step 7: Create quadrant classification
quadrant <- vector(mode = "numeric", length = nrow(localMI))
quadrant[DV < 0 & LM_I > 0] <- 1   # Low-Low
quadrant[DV > 0 & LM_I < 0] <- 2   # High-Low
quadrant[DV < 0 & LM_I < 0] <- 3   # Low-High
quadrant[DV > 0 & LM_I > 0] <- 4   # High-High
quadrant[localMI[, 5] > signif] <- 0  # Not significant

# Step 8: Assign quadrant to data
world_2024$quadrant <- quadrant

12.3 Plotting the LISA Cluster Map

The LISA Cluster Map provides a powerful visualization of localized spatial autocorrelation in happiness scores for 2024. By computing the local Moran’s I for each country, we can identify areas where similar values cluster together (e.g., high-scoring countries surrounded by other high scorers), or where outliers emerge (e.g., a low-scoring country surrounded by high scorers).

Each country is classified into one of five categories based on the direction and significance of the spatial relationship:

Category Description Statistically significant? Interpretation
HH High with high neighbours Yes Cluster (hot spot)
LL Low with low neighbours Yes Cluster (cold spot)
HL High with low neighbours Yes Spatial outlier
LH Low with high neighbours Yes Spatial outlier
Insignificant No strong autocorrelation No Insignifcant
Code 
# Load libraries
library(tidyverse); library(sf); library(spdep); library(tmap); library(readr); library(rnaturalearth)

# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
  mutate(year = as.numeric(year)) %>%
  filter(year == 2024 & !is.na(ladder_score))

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

# Join spatial and happiness data
world_2024 <- world %>%
  left_join(happiness, by = c("name" = "country")) %>%
  filter(!is.na(ladder_score))

# Construct spatial weights
nb <- poly2nb(world_2024, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Compute local Moran's I
localMI <- localmoran(world_2024$ladder_score, lw, zero.policy = TRUE)

# Prepare LISA categories
quadrant <- vector(mode = "numeric", length = nrow(localMI))
lagged_score <- lag.listw(lw, world_2024$ladder_score)
centered_lag <- lagged_score - mean(lagged_score)
centered_localMI <- localMI[, 1] - mean(localMI[, 1])
significance_level <- 0.05

# Define cluster categories
quadrant[centered_lag < 0 & centered_localMI > 0] <- 1  # Low-Low
quadrant[centered_lag > 0 & centered_localMI < 0] <- 2  # Low-High
quadrant[centered_lag < 0 & centered_localMI < 0] <- 3  # High-Low
quadrant[centered_lag > 0 & centered_localMI > 0] <- 4  # High-High
quadrant[localMI[, 5] > significance_level] <- 0        # Not significant

# Add to spatial data
world_2024$quadrant <- quadrant
world_2024$cluster_label <- factor(
  quadrant,
  levels = 0:4,
  labels = c("Insignificant", "Low-Low", "Low-High", "High-Low", "High-High")
)

# Plot the interactive LISA cluster map using tmap
tmap_mode("view")

tm_shape(world_2024) +
  tm_fill(
    col = "cluster_label",
    palette = c(
      "Insignificant" = "#ffffcc",
      "Low-Low" = "blue",
      "Low-High" = "#78c679",
      "High-Low" = "#c2e699",
      "High-High" = "red"
    ),
    title = "LISA Cluster (2024)",
    style = "cat",
    id = "admin",  # <<-- This makes popups/hover show country names
    popup.vars = c(
      "Country" = "admin",
      "Cluster Type" = "cluster_label",
      "Happiness Score" = "ladder_score"
    )
  ) +
  tm_borders(alpha = 0.4) +
  tm_layout(
    frame = FALSE,
    legend.outside = TRUE
  )

12.4 User interface design - Combining Moran Scattterplot with LISA map

Code 
# Load libraries
library(shiny); library(tidyverse); library(sf); library(spdep); library(tmap); library(readr); library(rnaturalearth)

# Load and prepare data
happiness <- read_csv("data/world_happiness.csv") %>%
  mutate(year = as.numeric(year)) %>%
  filter(!is.na(ladder_score))

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

# UI
ui <- navbarPage("World Happiness Explorer",

  tabPanel("Spatial Clustering",
    fluidPage(
      titlePanel("🌐 Spatial Clustering Analysis"),

      sidebarLayout(
        sidebarPanel(
          selectInput("selected_year", "Select Year:",
                      choices = sort(unique(happiness$year)),
                      selected = 2024),
          hr(),
          h4("Chart Interpretation"),
          helpText("The Moran scatterplot shows how each country's happiness score correlates with its neighbors'."),
          helpText("The LISA Cluster map highlights statistically significant spatial clusters:",
                   "\n- High-High: Happy countries near other happy countries",
                   "\n- Low-Low: Unhappy countries near unhappy neighbors",
                   "\n- High-Low / Low-High: Potential outliers",
                   "\n- Insignificant: No strong spatial pattern")
        ),

        mainPanel(
          fluidRow(
            column(6, plotOutput("moran_plot", height = "500px")),
            column(6, tmapOutput("lisa_map", height = "500px"))
          )
        )
      )
    )
  )
)

# Server
server <- function(input, output, session) {

  world_data <- reactive({
    data <- happiness %>%
      filter(year == input$selected_year)

    world %>%
      left_join(data, by = c("name" = "country")) %>%
      filter(!is.na(ladder_score))
  })

  # Spatial weights and local moran
  local_moran <- reactive({
    data <- world_data()
    nb <- poly2nb(data, queen = TRUE)
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
    localMI <- localmoran(data$ladder_score, lw, zero.policy = TRUE)

    lagged_score <- lag.listw(lw, data$ladder_score)
    centered_lag <- lagged_score - mean(lagged_score)
    centered_localMI <- localMI[, 1] - mean(localMI[, 1])

    quadrant <- vector(mode = "numeric", length = nrow(data))
    quadrant[centered_lag < 0 & centered_localMI > 0] <- 1  # Low-Low
    quadrant[centered_lag > 0 & centered_localMI < 0] <- 2  # Low-High
    quadrant[centered_lag < 0 & centered_localMI < 0] <- 3  # High-Low
    quadrant[centered_lag > 0 & centered_localMI > 0] <- 4  # High-High
    quadrant[localMI[, 5] > 0.05] <- 0                      # Not significant

    data$quadrant <- quadrant
    data$cluster_label <- factor(
      quadrant,
      levels = 0:4,
      labels = c("Insignificant", "Low-Low", "Low-High", "High-Low", "High-High")
    )

    list(data = data, lw = lw)
  })

  # Moran plot
  output$moran_plot <- renderPlot({
    dat <- local_moran()
    moran.plot(dat$data$ladder_score, dat$lw,
               labels = dat$data$name,
               xlab = paste("Happiness Score (", input$selected_year, ")", sep = ""),
               ylab = "Spatially Lagged Happiness Score",
               zero.policy = TRUE)
  })

  # LISA map
  output$lisa_map <- renderTmap({
    tmap_mode("view")

    tm_shape(local_moran()$data) +
      tm_fill(
        col = "cluster_label",
        palette = c(
          "Insignificant" = "#ffffcc",
          "Low-Low" = "blue",
          "Low-High" = "#78c679",
          "High-Low" = "#c2e699",
          "High-High" = "red"
        ),
        title = paste("LISA Cluster (", input$selected_year, ")", sep = ""),
        style = "cat",
        id = "admin",
        popup.vars = c(
          "Country" = "admin",
          "Cluster Type" = "cluster_label",
          "Happiness Score" = "ladder_score"
        )
      ) +
      tm_borders(alpha = 0.4) +
      tm_layout(frame = FALSE, legend.outside = TRUE)
  })
}

shinyApp(ui = ui, server = server)

Shiny applications not supported in static R Markdown documents